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The modulational instability and discrete matter wave solitons in dipolar BEC, loaded into a deep 
optical lattice, are investigated analytically and numerically. The process of modulational instability 
of nonlinear plane matter waves in a dipolar nonlinear lattice is studied and the regions of instability 
are established. The existence and stability of bulk discrete solitons are analyzed analytically and 
confirmed by numerical simulations. In a marked contrast with the usual DNLS behavior (no dipolar 
interactions), we found a region where the two fundamental modes are simultaneously unstable 
allowing enhanced mobility across the lattice for large norm values. To study the existence and 
properties of surface discrete solitons, an analysis of the dimer configuration is performed. The 
properties of symmetric and antisymmetric modes including the stability diagrams and bifurcations 
are investigated in closed form. For the case of a bulk medium, properties of fundamental on-site and 
inter-site localized modes are analyzed. On-site and inter-site surface localized modes are studied 
finding that they do not exist when nonlocal interactions predominate with respect to local ones. 
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I. INTRODUCTION 

The dynamics of Bose-Einstein condensates (BECs) in 
optical lattices has been the subject of intensive theoreti- 
cal and experimental investigations in recent times [TJ [2] . 
Different phenomena like generation of coherent packets 
of matter waves (atom lasers) [3], Bloch oscillations [IHB], 
gap solitons [7], discrete breathers [8, 9 j, compactons [10] 
have been predicted and experimentally observed. 

The possibility to vary different parameters of BEC 
systems - trap potential and atomic interactions - makes 
BEC a unique system for modeling different fundamental 
phenomena in condensed matter physics. The control of 
the strength and shape of the trapping potential induced 
by counter-propagating laser fields, as well as the time 
modulations of the parameters of the optical lattice, can 
be easily achieved. The strength of atomic interactions 
and, thus, the mean field nonlinearity can also be tuned 
in space or time by using, for example, the Feshbach res- 
onances method [11]. This allows the atomic scattering 
length a s to be tuned by a time-or-space variation of the 
external magnetic field near the resonant value. 

Joint effects of nonlinearity, periodicity and quantum 
pressure, lead to the existence of stable localized states 
conserving the form upon propagation and collisions. 
Gap solitons, discrete breathers and discrete compactons 
are examples of such states. Some of these structures are 
observed in experiments. Nonlinearity has usually been 
considered as local. Recent discovery of dipolar conden- 
sate with long-range interactions between atoms, posed 
the question of the existence of solitons in a dipolar BEC 
loaded in the optical lattice. The existence of solitons 
is due to the interplay between local (contact interac- 
tions) and nonlocal (long range interactions) nonlineari- 



ties and the optical lattice effects. Previously, the exis- 
tence of discrete solitons in systems with nonlocal interac- 
tions has been studied for semiconductors amplifiers [12] 
and nematic liquid crystals [13]. The first observation 
of a condensate in a gas of chromium atoms ( 53 Cr) with 
long-range interactions was reported in Ref. [14] . Dipolar 
condensate exhibits many unusual properties not encoun- 
tered in BECs with local interactions, e.g. the existence 
of stable isotropic and anisotropic 2D solitons [T5] , 

Two limiting cases can be distinguished: shallow and 
deep optical lattices. Both cases have been considered re- 
cently [16]. In the case of deep lattices, a discrete model 
using a nonlocal Gross-Pitaevskii equation is investigated 
in Ref. [T7] , Here, the dynamics of unstaggered bright 
discrete solitons for a 2D dzs£;-shaped dipolar BEC was 
analyzed. The system was studied by a 2D model based 
on the discrete nonlinear Schrodinger (DNLS) equation, 
with on-site and long-range cubic nonlinearities. The ex- 
istence of stable fundamental discrete solitons of the dif- 
ferent symmetries was shown. The authors also observed 
a stability exchange of the fundamental solutions but did 
not study in detail the region where both solutions are si- 
multaneously unstable. This is one the main goals of the 
present work, to explore in detail these regions and their 
dynamical properties. We will show that a very good mo- 
bility can be observed for solutions with a large value of 
the norm, what is absolutely forbidden in systems with 
only local nonlinearities. 

It is currently of interest to investigate discrete 
breathers in quasi-one-dimensional czgar-shaped dipolar 
BEC loaded in a deep optical lattice. Recently, this prob- 
lem was considered in Ref. [15] , where a non-polinomial 
DNLS model was obtained. The existence and stability 
of unstaggered bulk bright breathers was studied. Be- 
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sides unstaggered bulk solitons, there is especial interest 
on staggered bulk solitons, existence of surface discrete 
breathers (solitons) and the modulational instability of 
matter waves in the dipolar BEC embedded in an op- 
tical lattice. Surface solitons are the generalization of 
nonlinear Tamm states, and they have been well-studied 
in DNLS systems with on-site nonlinearity [TQl [20] . The 
existence of nonlinear Tamm states in quantum dipolar 
gases in optical lattices is a fundamental problem. We 
can expect a rich variety of Tamm states due to the com- 
petition between local on site and nonlocal cubic interac- 
tions and the lattice potential. In the present work, we 
explore in detail this issue. 

This paper is organized as follows: Section II intro- 
duces the discrete model of a Bose-Einstein condensate in 
a deep optical lattice subjected to dipolar interactions; in 
section III, we examine the modulational stability prop- 
erties of nonlinear plane matter wave solutions; in section 
IV, we review some results for bulk localized modes, in- 
cluding some newly found bistable behavior that contrast 
markedly with bulk phenomenology found in usual DNLS 
systems. In section V, we focus on surface localized states 
and, finally, section VI concludes the paper. 

II. LATTICE MODEL 

The system under consideration is the quasi one- 
dimensional dipolar BEC in a deep optical lattice. The 
governing equation is the ID Gross-Pitaevskii equation 
(GPE) with nonlocal interaction terms [21]: 

d^r h 2 d 2 ^f 

+ *(x,t) / d£R(\X-£\) m,t)\ 2 = 0, (1) 

where gin = 2ha s uj_. u±_ corr esponds to the transverse 
trap frequency, l± = ^/h/muj_\_, and d is the dipolar mo- 
ment. The parameter a can vary from 1 to —1/2 for 
dipoles oriented along or perpendicular to the x-axis. 
The wave function is normalized to the number of atoms 
comprising the BEC, N = \^(x)\ 2 dx. Now we define 
dimensionless parameters: 

h 2 k 2 V 
E R = — — = hbj R = — , t = Tuj r , x = kX, 
2m V 

aa d a s l2a s uj± 

9o = i—j , Qo = — , ^ = \ ®, 

l_Lka s0 a s0 V ur 

where = md 2 d /h 2 is the characteristic scale of the long- 
range dipolar interactions, and a s o is the background 
value of an atomic scattering length. Equation ([I]) is 
recasted as: 

/OO 
dZR(\x-Z\) \^,t)\ 2 = (2) 
-OO 



where ip(x,t) represents the mean-field wave function of 
the condensate. Two commonly used nonlocal kernel are, 

R x {x) = (l + 2x 2 )exp(x 2 )erfc(|x|) -2tt" 1/2 |x|, and 

R 2 (x) = 8\x 2 + 8 2 )-^ 2 . 

The first kernel corresponds to a dipolar BEC in a quasi- 
1D trap [21 ; while the second kernel, with the cutoff 
parameter 5, is more convenient for an analytical treat- 
ment [22 . The matching condition i?i(0) = i?2(0) im- 
plies S = 7r -1 / 2 . Parameter S corresponds to the effective 
size of the dipole. It takes a value of the order of the 
transverse confinement length, which makes the model 
one-dimensional, and constitutes the unit of length in 
Eq.|l]). Therefore, the choice of 0.56 is quite 

reasonable. In the limit x ^> 5, where dipole-dipole in- 
teraction effects dominate the contact interaction effects, 
both response functions behave as ~ 1/x 3 . 

In a deep optical lattice, V ^> 1 and it is natural to 
consider the expansion 

OO 

n= — oc 

where <j) n (x) are Wannier-like functions located on the 
minima of the periodic potential V(x). The analysis of 
overlap integrals [23, 24 shows that the equations for u n 
become: 

g(\u n+ i\ 2 + |^ n _i| 2 )^ n = 0, (3) 

where 

Q = Qo J J R(x - i)\(j) n (x)\ 2 \(j) n (i)\ 2 dx d^ and 

9 = 90 J J R(x-0\^n±i(x)\ 2 \^ n (C)\ 2 dxdC • 

It should be noticed that the parameter q, which is pro- 
portional to the atomic scattering length a 3 (t), can be set 
to zero by means of the Feshbach resonance method [11] . 
According to this technique, by a variation of the exter- 
nal magnetic field near the resonant value, it is possible 
to diminish to zero the atomic scattering length respon- 
sible for the mean field nonlinearity parameter q. This 
limit of the model has been applied recently to the analy- 
sis of the fundamental limit on the atomic interferometer 
based on BEC with tunable scattering length, loaded in 
optical lattices [24]. 

Equation (|3| possesses two conserved quantities, the 
norm N and the Hamiltonian H, 

M 

n = ^Kl 2 

n=l 

M 

H = [^n+l< + ||^n| 4 + ||^n+l| 2 |^n| 2 +C.C. 

n=l 
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where M indicates the number of lattice sites. Linear 
plane- wave solutions of Eq.(3| take the form u n (t) = 
uq ex.-p(ikn + ujt) and satisfy the linear dispersion relation 
uo — 2k cos fc, which defines the band of single-particle en- 
ergies uok G {—2ft, 2ft}. Outside of this band, nonlinear 
localized solutions are expected to exist, uo and k corre- 
spond to the normalized amplitude and quasimomentum 
of the condensate, respectively. 

Along this work, we will look for bulk and surface one- 
dimensional fundamental, centered on-site and inter-site, 
solutions. By implementing a standard Newton- Raphson 
method, we numerically compute localized stationary so- 
lutions of model (|3|, of the form u n (t) = u n exp [iut], by 
solving the set of equations 

uo u n = (u n +i + u n _i) + qui + g(u 2 n+1 + *4-i)^n- (4) 

where u n G R. Hereafter we will consider ft = 1 and 
q > 0. The regime q < can be explored by simply 
making the transformation u n — ^ (— l) n u ni uo — » —uo and 
g — >• —g. For each solution, we characterize it by comput- 
ing its norm and Hamiltonian. We perform a linear sta- 
bility analysis by using a standard method developed in 
Ref . [25] . We obtain an eigenvalue spectrum and compute 
the eigenvalue with the largest imaginary part (denoted 
as G) as an indication of the instability gain. When plot- 
ting the norm- frequency diagrams of each mode (bulk or 
surface), we will use continuous (dashed) lines to denote 
stable (unstable) solutions. In addition, and as a visual 
aid, quantities of interest for the on-site (inter-site) solu- 
tions will be plotted in black (orange). 



leading to the linear system 

[—ft — uo(k) + a + ] u\ J rbu2 
b u\ + [O — uo(k) + a~) U2 



(8) 



where b = 

2cos(Q±/c) 

vided 



qui 



2guQCOsQ and a ± = 2(q + g)ul 



2gul cos Q. Nontrivial solution exists pro- 



d± v / d 2 - 4b 2 + 4(uj(k) - a+)(uo(k) - a~) 



(9) 

where d = a + — a . From this, we define the instability 
gain G = Im[fi]. If G ^ 0, the nonlinear plane wave 
experiences "modulationally instability" (MI). 

Let us discuss the stability of uniform plane waves 
(k = 0). In the special case g = 0, the system re- 
duces to a pure DNLS chain [27] and the gain for the 
plane wave becomes G = ±Im[y^2 sin(Q/2) 2 — quQ, as 
shown in FigfTJa), where G is shown in the form of a 
density plot as a function of uo and Q. For < qu^ < 2, 
there is MI for Q < 2 arcsin( ^/qul/2). For qu^ > 2 
we have MI for all Q values. Therefore, plane matter 
waves will be always unstable. In the case of zero effec- 
tive mean field nonlinearity (q = 0) the gain is given by 
G = 4Im{[(cosQ-l)(-l + (l + 2^)cosQ)] 1 / 2 }. Anal- 
ysis of this expression shows that, for g > and a given 
uo, there is always a Q-interval where G / 0, and the 
uniform plane wave is always unstable. On the contrary, 
for g < 0, the plane wave is stable for \uo\ 2 < l/\g\. For 
q,g > 0, it can be proven from Eq.(|9]) that G =^ for 



III. MODULATIONAL STABILITY 

We begin by studying the linear stability of plane 
waves solutions under the effect of nonlinear interactions. 
As pointed out a long time ago [26], modulational insta- 
bility in discrete systems is a very efficient mechanism 
to generate discrete solitons. Equation ([3| admits plane- 
wave solutions that lead to the nonlinear dispersion rela- 
tion 



uo(k) = 2 cos k + (q + 2g)u 2 ) . 



(5) 



Next, we compute the modulation stability of this plane 
wave solution by setting a perturbed solution in the form 
u n (t) = [uq + Suo(t)] ex.p(ikn + cut). After replacing in 
Eq.Q and keeping linear terms in Su n only, we obtain 
the evolution equation for the perturbation, 

(5u n+1 e ik + * n _i exp- i/c ) + qi%5u* n + (6) 
gul(5un+i + Su n -i + + ^<_i) = . 



We pose Su n (t) in the form 

5u n (t) = ^e i(gn+m) 
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Uo 1 -0 
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tt/4 tt/2 3tt/4 tt tt/4 tt/2 3n/4 n 




* -i(Qn+Q*t) 



u 9 e 



(7) 



FIG. 1: (Color online) G for nonlinear plane wave solution 
for q = 1 and g = (a), g = 1 (b), g = -0.2 (c), g = -1 (d). 
Darker (lighter) color means a stable (unstable) region. 
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any uq and, therefore, the plane wave is always unstable 
[see Fig. [TJb)] . For q > 0,g < 0, it can be shown that 
for \g\ < q/2, G / 0, while for \g\ > q/2, G = for 
ul < 2/(q + 2|#|) [see Figs. 0c) and (d)]. 

The above results suggest that formation of discrete 
solitons will be likely in most cases, with the exception of 
strong attractive dipolar interactions, where a minimum 
norm will be required. These rough predictions will be 
confirmed in Section V. 



IV. DIMER APPROACH 

In order to get a deeper understanding of the present 
model, it will prove useful to consider first the dimer limit 
M = 2 [23] . This constitutes an integrable system which 
can give us some insights about the general phenomenol- 
ogy occurring in larger systems. In particular, it will 
prove useful when we consider localized surface modes. 
We solve Eq.Q for M = 2 



UJU\ 
UJU 2 



U 2 
Ul 



qu\ 



gu\u x 



qu\ + gu\u 2 



(10a) 
(10b) 



As a general ansatz, we consider u\ — A and U2 = ft A. 
After inserting this ansatz in (10), we obtain three sta- 
tionary solutions: 



P = ±1 and (3 = 



1 



(q-g)A* 



The solution /3 = 1 corresponds to a "symmetric" so- 
lution (ui = U2) for which N sym = 2(u — l)/(q + g). 
This solution bifurcates from the symmetric linear mode 
at (j = 1. A second solution, /3 = — 1, corresponds to 
the "antisymmetric" mode {u\ = —U2) with N ant = 
2(u + l)/(q + g), bifurcating from the antisymmetric lin- 
ear mode at uj = — 1. A third solution is called "asym- 
metric" because, in general, \u\\ ^ \u2\. It appears at 
oJmin = 2q/\q - g\, the only frequency where \m\ = \u 2 \, 
bifurcating from the "symmetric" ("antisymmetric") so- 
lution if q > g (q < g) [see thick orange lines in Figj2]. 
This also implies a change in the solution's topology from 
"in-phase (ip)" to "out of phase (op)". For this asym- 
metric mode N asy = u/q, i.e., this mode does not depend 
at all on the dipolar interaction. A standard linear sta- 
bility analysis show that the symmetric solution is stable 
for N < 2/(q — g) [i.e., before the onset of the stable 
asymmetric solution, for q > g] while the antisymmetric 
mode is stable for N > 2/(g — q) [i.e., after the onset of 
the unstable asymmetric solution, for q < g]. By fixing q 
while increasing g, we observe that the slope of the N vs. 
uj curves for the symmetric and antisymmetric solutions 
decreases (for q — g both curves coincide with the asym- 
metric one). It is important to notice that the symmetric 
solution increases its stable existence region while uJmin 
increases as g —> q. For g > q, the symmetric solution 
is always stable; the asymmetric one becomes unstable 
bifurcating, now, from an always stable antisymmetric 
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FIG. 2: (Color online) Norm versus frequency diagrams for 
q — 1. Black and gray thin lines correspond to the symmetric 
and antisymmetric solutions, respectively, (a) g = 0.5 (con- 
tinuous) and g — 1.5 (dashed), (b) g — —0.5 (continuous) and 
g — —1.5 (dashed). The thick line represents the asymmetric 
solution. Insets depict stationary dimer profiles. 



solution (for q > 0). Figj2|a) shows an example of this 
phenomenology, where uj m i n = 4 for g = 0.5 and g = 1.5. 
For the sake of simplicity, we have plotted the asymmet- 
ric solution as a "continuous" thick orange line for both, 
the stable and unstable, cases. 

We now decrease g from zero while keeping q fixed. We 
see that when g —> —q, the norm of the symmetric and 
antisymmetric solutions diverges. Before this happens, 
the global sign of the nonlinearity (q + g > 0) and the 
slope is positive. However, when q + g < the effective 
sign becomes negative and the slope changes completely. 
The asymmetric solution increases its norm as before, be- 
cause its slope does not depend at all on the g- value; i.e., 
as soon as q > this solution possesses a positive slope 
with uJmin > for any g. See Fig{2|b) as an example of 
this phenomenology for g = —0.5 and g = —1.5. From 
([5| we can see that, for a larger system, this situation 
will occur when g — » —q/2, due to the larger number of 
nearest-neighbors. For g < —q/2, the fundamental lin- 
ear mode (k = 0) will increase its norm by decreasing 
its frequency. Of course, this will have consequences for 
localized solutions bifurcating from the linear band, as 
we will see below. 



A. Effective potential 

We now construct an effective potential for the dimer 
model that connects the stationary solutions found above 
with a dynamical picture of this problem [32]. This will 
be important to better understand the method we will 
implement below when dealing with larger systems. For 
any stationary solution, we can define a center of mass 



as p = u\j~N being N = u\ 



x 2 . 



p = and p = 1 



denote a solution located at site n = 1 and site n = 2, 
respectively; while p = 0.5 implies a solution centered at 
the inter-site position (|/3| = 1). We now express u\ and 
U2 in terms of the p coordinate: u\ = ±^iV(l — p) and 
U2 = ±y/Np). After inserting these expressions in the 
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FIG. 3: Effective potential for q = 1, g — 0.5 (left) and g — 1.5 
(right); for N — 1 (black line) and N = 10 (gray line), if has 
been scaled for comparison purposes. 



dimer Hamiltonian, we obtain 



H(p) = ±2N^p(l^p) — N 2 ± + (g- q)p(l - p) 



where — (+) denotes the ip (op) cases. For a given norm 
N, we look for critical points of the effective potential 
H(p), obtaining 



Psym,ant 



and 



Pasy — ^ ^ 



1 



as critical points. Therefore, if TV > 2/\q — g\, there 
are always three stationary solutions: one centered in 
between the two sites (symmetric or antisymmetric) and 
two asymmetric solutions with a varying center of mass. 
For N min = 2/\q — g\, the asymmetric solution simply 
coincides with the symmetric or antisymmetric one [black 
dots in Figj3] where N min = 4]. As the norm increases, 
one asymmetric solution bifurcates to the left and the 
other to the right from p = 0.5 [see gray points in Figj3]. 
For N ^> N min , the asymmetric solutions locate at p — > 
0,1. 

The sketched potentials agree perfectly with the prop- 
erties described before. For g < q and N < N mini the 
only stationary solution is the symmetric one which cor- 
responds to a minima (stable) in the effective potential 
[see black line in Fig j3]- left] . Asymmetric solutions exist 
only above some norm threshold and they correspond to 
two local minima (stable) in the effective potential which, 
as expected, possesses a local maxima in between corre- 
sponding to the unstable symmetric solution. Therefore, 
the appearance of the asymmetric solution destabilizes 
the symmetric one. On the other hand, for g > q the 
situation is quite different. For a small norm the anti- 
symmetric solution is a maximum (unstable) of the effec- 
tive potential [see black line in Figj3]-right] . For a larger 
norm, the asymmetric solution appears as a maximum in 
the potential H(p) being, therefore, unstable while the 
antisymmetric solution stabilizes. This means that, in 
this case, the appearance of the asymmetric solution sta- 
bilizes the antisymmetric one. 



V. BULK LOCALIZED SOLUTIONS 

In this section, we will study localized stationary so- 
lutions located at the center /middle of the lattice (bulk 
solitons). Some profiles are shown in the inset of Fig{4] 
for a lattice of M = 100 sites. In general, the profiles 
change in a very defined manner. For a fixed norm we 
see that, by increasing the value of g from zero, on-site 
profiles increase its width by increasing the amplitude 
of the first nearest-neighbors while the inter-site modes 
does not change significantly. Therefore, a positive dipo- 
lar long-range interaction promotes a wider on-site profile 
due to the increasing relevance of nonlocal amplitudes. 
For the same reason, the inter-site modes are not effec- 
tively affected by this increasing g. On the other hand, 
for a negative value of g, solutions become more local- 
ized. A decreasing long-range interaction makes the first 
nearest-neighbor amplitudes smaller, leading to profiles 
that are very well localized in space. 



A. q,g>0 

First of all, we consider the case where both nonlin- 
ear coefficients are positive. We also fix, without loss of 
generality, q = 1 and vary g > 0. For g = 0, we obtain 
the well-known DNLS phenomenology where the on-site 
solution is always stable while the inter-site mode is al- 
ways unstable [thick lines in Figj4j . If we increase the 
value of g, for instance to g = 0.8 [thin lines in FigJI], 
we see how the slope of the inter-site family decreases 
abruptly, getting closer to the one for the on-site solu- 
tion. (This is somehow similar to the situation encoun- 
tered for the dimer, where the slope of the symmetric 
solution approaches the one of the asymmetric solution). 
In addition, the stability properties start to change for g 
approaching q. At low frequencies, there is a small re- 
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FIG. 4: (Color online) Norm versus frequency diagrams for 
q = 1. Black (orange) curve represents the on-site (inter-site) 
mode. Thick (thin) lines correspond to g = (g = 0.8). Inset: 
Profiles at N = 25 for g = (thick) and 0.8 (thin). 
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gion where the on-site solution becomes unstable, while 
the inter-site mode stabilizes. Then, both solutions ex- 
change stability and the previously described behavior is 
reversed. However, if we increase the value of g, a bit 




FIG. 5: (Color online) Norm versus frequency diagram for q — 
1 and g = 0.96. Black (orange) thick lines correspond to on- 
site (inter-site) modes, while black thin lines to intermediate 
solutions. Inset: Instability gain versus frequency. 

more, for instance to g = 0.96 [see Fig{5], we first see 
- in the N versus uj diagram - how the solution fami- 
lies approach, being almost undistinguishable from each 
other. Now, if we take a look of the stability properties, 
a very interesting phenomenology emerges. The previ- 
ous apparently trivial exchange of stability is not such. 
For lower frequencies, as described before, the inter-site 
solution is stable while the on-site one is unstable; i.e, 
a completely opposite stability picture compared to the 
DNLS limit. Then, there is a complete region where both 
fundamental solutions are simultaneously "unstable" [see 
Figj5]-inset] . For discrete nonlinear systems this kind of 
behavior was also observed in more complex models [28]- 
[30] which include the same linear and nonlinear disper- 
sion terms plus many others. On the other hand, a com- 
pletely opposite phenomenology have been predicted for 
optical saturable systems [3TJ [32] with multiple regions 
of simultaneously "stable" solutions. Therefore, from the 
fundamental and dynamical point of view, the present 
phenomenology is very important and, certainly, inter- 
esting because, in principle, we could expect good condi- 
tions for mobility in one and, also, higher dimensions. 

In order to delve deeper in this analysis, we study in 
detail the stability of both fundamental solutions [see 
Figs(6ja) and (b) where a brighter (darker) color cor- 
responds to a more unstable (stable) solution]. First of 
all, from this analysis it is clear that the on-site (inter- 
site) solution is stable (unstable) if g < q, while for g > q 
the situation is the opposite. Figj6|c) shows the region, 
in parameter-space, where both solutions are "simultane- 
ously" unstable. As g increases the exchange region also 
increases showing that, by tuning the dipolar interaction, 
we could observe this phenomenology for different values 



of the norm. We computed the width (A7V) of this bi- 
unstable region [see FigjHJd)] for different values of g. We 
can see that A TV tends to diverge when g — >• indicating 
the full stability exchange for g > q. From a dynamical 
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FIG. 6: (Color online) Instability gain as a function of g and 
N for (a) on-site and (b) inter-site solutions, (c) Bi- unstable 
region, (d) g versus AiV for the bi-unstable region, (q = 1). 

point of view, theory tells us that once we have two un- 
stable solutions, both should correspond to a local max- 
ima in a Hamiltonian representation. Therefore, there 
should exist another solution in between corresponding 
to a local minima, a stable "intermediate" solution (IS) 
[see inset in Figj7]. The IS is an asymmetric station- 
ary solution with a profile that varies from an inter-site 
mode (smaller norm) to an on-site mode (larger norm). 
In that sense, the IS is a kind of stability carrier that 
exchanges the stability of both fundamental solutions. 
Now, we compute an effective potential, i.e. the Hamil- 



tonian (H) versus the center of mass (p = ^ n \ 



i N ) 



of different solutions across the lattice, by using a con- 
straint method [20] [3Tj similar to the computation per- 
formed for the dimer. (An on-site solution possesses an 
integer p-value while the inter-site mode a semi-integer 
one). In this picture, a stable solution will correspond to 
a local minima while an unstable solution coincide with 
a local maxima. For example, in a DNLS phenomenol- 
ogy the on-site solution corresponds to a minima, while 
the inter-site one to a maxima [see the gray thin curve 
in Figj7|as an example of this behavior]. Therefore, the 
effective potential in cubic systems is periodic, and the 
energy differences between these two fundamental solu- 
tions constitute the energy barrier in the system (also 
known as the Peirls-Nabarro barrier). Mobility will only 
be possible if the solution has enough kinetic energy to 
overcome this self- induced energy barrier [33]. However, 
in more complex systems, the situation is not that simple. 
For instance, in the present model, the effective energy 
barriers will depend on the particular value of the norm, 
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g, and the particular stability region. FigjT] shows dif- 
ferent effective potentials by fixing the value of g = 0.96 
and varying iV, where we have plotted "normalized" H- 
values in order to compare the different cases. First, we 
can see that for a smaller norm (black line), the inter- 
site solution corresponds to a minima (stable) while the 
on-site mode corresponds to a maxima (unstable). In 
the region where both fundamental solutions are simul- 
taneously unstable (thick orange line), both correspond 
to a maxima. The stable IS corresponds to a minima 
located in between both fundamental solutions [see inset 
in Fig{7]. The potential is again periodic but, now, it has 
a more complicated "binary" geometry. The potential is 
softer when going from the IS to the inter-site one than 
when going to the on-site mode. The shape of this po- 
tential is also a consequence of the stability properties for 
both fundamental solutions. For this value of the norm, 
the on-site solution is more unstable than the inter-site 
one [see Fig J5} inset], what coincides with the shape of 
the two maxima. From this picture, we could predict 
that an on-site solution will require a smaller kinetic en- 
ergy to move than the inter-site mode. By increasing 
the norm, the stability analysis predicts that the on-site 
(inter-site) mode becomes stable (unstable). The effec- 
tive potential becomes a simple cubic-like potential as in 
the DNLS limit. 
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FIG. 8: (Color online) (a) and (b) Evolution of \u n \ 2 of an 
initial IS (N = 48) for k = and for k = 0.18, respectively, 
for tmax = 100. (c) Center of mass versus time for k — 0.18. 
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FIG. 7: (Color online) Hamiltonian versus center of mass por 
N = 35 (black line), N = 48 (orange thick line) and N = 59 
(gray thin line) for q — 1 and g — 0.96. Inset: Intermediate 
solution for N = 48. 



In order to corroborate the potential's picture, we in- 
tegrate numerically model (|3| to study the dynamics for 
N = 48. We initialize our numerical integration with 
the three stationary solutions of this problem includ- 
ing a kinetic energy with a 'kicked' mode of the form: 
^n(0) = u n exp [ik(n — n c )]. First of all, we observe a 
stable propagation of the intermediate solution (IS) [see 
Figj8^a)], as expected from the computed potential and 
the stability analysis. Contrary to what is typically ex- 
pected, this asymmetric solution propagates stably for 
long times. From Fig(7|we know that some energy should 
be added to the IS in order to set it in motion. Fig ji 



shows an example of very good mobility with no vis- 
ible radiation from tails; i.e., a coherent movement of 
the atomic wave-packet. The solution just propagates 
"feeling" the topology of its own potential. Figj8^c) 
shows how the velocity (slope) of the propagating so- 
lution changes according to the energy surface. When 
it passes through integer (thin horizontal lines) or semi- 
integer (dashed horizontal lines) positions, the velocity 
decreases because these points correspond to local max- 
ima. On the other hand, when it passes through positions 
related to IS (dotted horizontal lines), the velocity in- 
creases, as expected from general dynamical arguments. 

Now, we initialize our numerical simulation with un- 
stable on-site and inter-site solutions at N = 48. It is im- 
portant to notice that, in this case, H onsite > H intersite 
[see FigjT]. The on-site solution is located in a very 
sharp local maximum, so a very unstable dynamics is 
expected. In addition, as this maximum is larger than 
the one for the inter-site solution, a good mobility is ex- 
pected through the lattice, if radiation from tails is not 
too large. Figs{9ja) and (c) show the propagation of an 
on-site solution without any initial kinetic energy. Dy- 
namics corroborate the expected phenomenology from 
the potential analysis: the velocity tends to zero when 
the solution passes through an on-site position (p = 51, 
52 and 53, in this example); the velocity is small but not 
zero when solution crosses an inter-site position [p = 51.5 
and 52.5, in this example); and the velocity increases to 
a maximum when the solution passes through an IS po- 
sition (p = 51.2, 51.8, 52.2 and 52.8, in this example). 



8 




50 



52 



54 



56 




53.0 
52.5 
52.0 
51.5 
51.0 



51.0 
50.9 
50.8 
50.7 
50.6 
50.5 











d) 











FIG. 9: (Color online) (a) and (b) Evolution of \u n \ for an 
initial on-site and inter-site solution, respectively, for N — 48. 
(c) and (d) Center of mass versus time for cases (a) and (b). 
k = and t m ax — 100. 



On the other hand, by initializing the numerics with an 
unstable inter-site solution, without any initial kinetic 
energy, the dynamics is less unstable than before. The 
solution takes more time to destabilize and to start an os- 
cillatory dynamics inside the potential well [see Figsj9^b) 
and (d)]. The profile changes from an inter-site mode 
(zero velocity) to an IS (maximum velocity) and then 
it approaches the on-site solution (zero velocity). Then, 
the solution goes back and the cycle starts again. This 
system could be used as an "atomic clock", where the 
fundamental period would depend only on the amount of 
particles in the system (N). 

For g > q, the inter-site modes become always stable 
while the on-site solutions are just unstable. This phe- 
nomenology is completely the opposite to the DNLS one, 
therefore when the dipolar interaction is larger than the 
mean field nonlinearity broader solutions are favored. 



B. q > 0, g < 

Now, we consider q = 1 and g < for which the inter- 
site mode is always unstable in the whole range of pa- 
rameters. For g < —0.3, the on-site solution presents 
a critical norm, where the slope changes its curvature 
and the solution destabilizes [see thick lines in FigflO^a)]. 
Then, it fuses with the inter-site mode becoming stable 
when changing again its curvature (similar to the case 
of a 2D DNLS model). After this, both modes just de- 
crease their norm up to the edge of the linear band at 
uj = 2. By decreasing the value of g further, we start to 
observe that the inter-site solution also shows a change 
of curvature for a given norm [see as an example thin 
lines in FigflO[a)]. Moreover, both fundamental solu- 
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FIG. 10: (Color online) Norm versus frequency diagrams for 
on-site (black) and inter-site (orange) modes, (a) g — —0.3 
(thick) and g = -0.48 (thin), (b) g = -0.5 (thick) and 
g = -0.58 (thin), (q = 1). 



tions still bifurcate from the linear band (uj — >• 2) with a 
large slope as g approaches the critical value of g = — 0.5. 
When g = —0.5 and, therefore, q + 2g = [see thick lines 
in FigflO^b)], both solutions increase their norm to in- 
finite at a given minimum norm value. The decreasing 
branch coming from uj ^> 2 should connect with a family 
which bifurcate from the linear band (uj ~ 2) at an in- 
finite rate, as predicted from Eq.([5|. Therefore, a norm 
threshold should appear in between in order to connect 
the two family branches. If we decrease the value of g 
even further [see thin lines in FigflO^b)], we observe that 
the solutions crosses the edge and enter into the linear 
band. There is a branch that originates at uj = 2 which 
increases with negative slope as Eq.([5]) predicts for the 
fundamental mode when q + 2g < 0. (This is similar to 
the behavior observed for the dimer in Figj2^b). In that 
case, the symmetric mode starts to increase its norm in 
the negative direction of frequencies, identical to the be- 
havior observed for bulk solutions originating from the 
linear band). For even more negative g- values, solutions 
start to deviate from each other with very different norm 
thresholds. The inter-site solution is always unstable and 
its slope and norm threshold increase as g decreases. In 
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fact, for g ^ —1 no inter-site solution was numerically 
found. On the other hand, the on-site solution does not 
feel the decreasing value of g, keeping its slope constant 
and being stable for all norms above norm threshold. 



VI. ANALYTICAL ESTIMATES 

In order to obtain a deeper understanding of our nu- 
merical findings, we test the expected phenomenology of 
the stationary model Q in two limit regimes: large and 
small norm. Let us consider first the large norm limit. 
In general, an on-site localized profile will always have 
a maximum amplitude A while its nearest-neighbor sites 
will have a value /3A, with |/3| < 1. On the other hand, an 
inter-site solution will always have two main peaks 
and two symmetric neighboring sites with amplitude eB 
(|e| < 1). By inserting these ansdtze in the equation for 
the center site, we get uj on-site — 2/3 + (q + 2gf3 2 )A 2 and 
Winter-site = (1 + e) + [q + g(l + € 2 )}B 2 . It is known 
that for large norm, solutions are, in general, highly lo- 
calized, i.e. /3, 6 —> 0. In addition, for an on-site mode 
the norm is essentially given by A 2 while for the inter-site 
mode is approximately 2B 2 . Therefore, to first approx- 

imation, N%? aite ~ u/q and A&- site ~ MGz + g). 
Thus, the increment in norm for a large frequency is lin- 
ear and independent of the g- value for the on-site mode. 
For the inter-site mode, there is an special case, q = g, 
where the two fundamental solutions coincide, with value 
N ~ oj/q. Therefore, for g < q, the inter-site solution 
has a larger norm for the same frequency, being there- 
fore unstable. However, for g > q the on-site solution 
has a larger norm for the same frequency, being now un- 
stable. In that sense, our estimate predicts a change in 
the stability properties for the fundamental solutions for 
q « g. On the other hand, for g < 0, our estimates pre- 
dict that the inter-site slope increases while the on-site 
one keeps equal. As an example, for g = —0.5 the inter- 
site norm is around four times larger than the on-site 
norm. If g — >• —1, our estimates say that the inter-site 
solution should diverge (disappearing for g < —1) while 
the on-site mode remain unaltered. All these analytical 
predictions are in perfect agreement with our numerical 
results shown in section IVl 

Let us take now the small norm limit. Typically, we 
know that when norm decreases solutions become wider, 
bifurcating from the linear band fundamental mode k = 
0, if q > (or k = tt if q < 0). Thus, in this limit, we can 
assume that /?, e — » 1 (meaning that all sites are approxi- 
mately equally excited) and N smal1 ~ M(u-2)/(q + 2g), 
for both fundamental solutions. This expression tells us 
that the initial slope is larger when the system is larger 
and that the initial frequency for this solution is uj = 2. 
These estimates also agree with our numerical findings, 
including the one which predicts that, for q-\-2g = 0, the 
solutions bifurcating from the linear band will diverge, 
generating a norm threshold that is independent of the 
lattice size. For q + 2g < our estimate suggests that 



uj should be smaller than 2 in order to keep the norm 
positive. Therefore, the solution bifurcating from the 
linear band increases its norm by initially decreasing its 
frequency. 

There is another issue related to a change of topology 
for the on-site solution. In the dimer problem, we saw 
that for g < q, the asymmetric solution was unstaggered 
(u\U2 > 0), same as the symmetric solution from where it 
bifurcates. However, for g > q, the asymmetric solution 
changed its topology becoming staggered (U1U2 < 0), 
same as the antisymmetric solution from where it bifur- 
cates. Let us try to predict what would happen for bulk 
solutions. We can do an analytic approximation of the 
profile by using "strongly localized modes (SLM)". We 
consider the following stationary ansatz for the on-site 
solution ~ A{0, ...,0,/3, l,/3,0, ...,0}. By replacing 
these expressions in Q, for the center and first nearest- 
neighbor sites, we find that 

R= -(q-g)A 2 ±^8 + (q-g) 2 A* 
p 4 

For q > g, the sign applies, being f3 > and the 
solution is unstaggered (u n u n +i > Vn). In addition, 
as g approaches f3 increases and the on-site solution 
becomes wider. For q < g the correct sign is "— ", imply- 
ing that /3 < and that the solution losses its topol- 
ogy. Therefore, similar to the dimer case, there is a 
change in the topology of solutions at least for q < g, 
where uj increases monotonically with the increment of 
the norm. We numerically found different on-site solu- 
tions for q < g. We continue the on-site unstaggered's 
family, however the first nearest-neighbor amplitudes be- 
comes very large and the SLM approach does not describe 
the solutions properly. On the other hand, we also found 
on-site solutions where the sign of the main peak was dif- 
ferent to the sign of the first nearest-neighbor amplitudes, 
coinciding with our prediction. However, this solution 
does not belong to the family of fundamental solutions 
we focus on this work. Finally, for g < (and q > 0), the 
on-site solutions are always unstaggered (/? > 0) being 
more and more localized as g decreases. 

For the inter-site solution we use the stationary ansatz 
u %n ^ ^ e ^ ^ e? o, 0}, obtaining 

_ -(l + g£ 2 ) + V4 + (l + ^ 2 ) 2 
6 " 2 

which displays no dependence on the g- value. That 
means that - for example - for q = this solution will 
not exist in the large norm limit (this is not the case of 
the on-site solution). Therefore, no topological transi- 
tions are expected for this type of solutions being always 
e > 0, for all g, g. 

VII. SURFACE LOCALIZED SOLUTIONS 

We consider now fundamental, on-site and inter-site, 
solutions located at the right surface of a ID lattice, i.e. 
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"surface solutions", as shown in Fig {TT]- inset. This kind 
of localized modes possesses a norm threshold for their 
excitation in usual DNLS lattices [20]. Therefore, it turns 
interesting to study the effect of the dipolar interaction 
in the excitation of such structures. 

The main phenomenology that we observe concerning 
norm thresholds, consists of their increment as g grows 
from zero. Therefore, it is more difficult to sustain a lo- 
calized solution at the surface if the long-range dipolar 
interactions increase. The increasing repulsive charac- 
ter of the surface become evident from the norm ver- 
sus frequency diagrams [see Fig 11 . In all these curves 
both - on-site and inter-site - solutions get connected af- 
ter norm threshold where the on-site solution changes its 
curvature and fuses with the inter-site mode. Therefore, 
when g — >> and the norm threshold tends to infinite, 
both solutions disappear. In this case, this happens when 
the inter-site solution decreases its slope and approaches 
the one for the on-site solutions, as Fig 11 shows. As 
a strong consequence, unstaggered on-site and inter-site 
surface solutions do not exist for g > q. 




FIG. 11: (Color online) Norm versus frequency diagrams for 
q — 1. Black (orange) curve represents the on-site (inter- 
site) modes. Thick lines corresponds to g = 0. Dashed and 
long-dashed thin lines correspond to g = 0.4 and g = 0.8, 
respectively. Inset: Profiles for N = 25 and for g = (thick) 
and 0.8 (thin). 

On the other hand, when g is lower than zero the norm 
threshold decreases. Fig 12 shows this behavior for g = 
—0.4 where both solutions connect each other as in the 
DNLS limit. However, for smaller g- values both solutions 
posses a norm/frequency threshold increasing their norm 
for lower frequencies. If g ^ —0.5, both solutions connect 
each other when approaching the linear band edge. For 
g = —0.5 both solutions tend to the linear band edge 
increasing their norm to infinite at uj — >• 2 [see Figjl2| . As 



2 [see Fig ; 

for the bulk solutions, the surface ones bifurcating from 
the linear band increase their norm because their minimal 
frequency decreases for g < —q/2 (the effective "linear- 
band-border" also decreases). Fig fl2] shows how, for g = 
—0.85, solutions touch the linear-band-border at uj = 2 
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FIG. 12: (Color online) Norm versus frequency diagrams for 
q — 1. Black (orange) curve represents the on-site (inter-site) 
modes. Thick lines corresponds to g — 0. Long, middle, and 
short dashed thin lines correspond to g = —0.4, g = —0.5 and 
g = —0.85, respectively. 



but at very different norm values. As we analytically 
predicted for inter-site bulk solutions, for surface ones 
the norm also diverges when g — >> —q. Therefore, the 
slope of this family increases very rapidly in comparison 
to the on-site one [see Fig 12 . That implies a larger 



norm threshold for the inter-site mode, which diverges 
when g — >• —q. Fig 12 shows that the large norm limit 
for the on-site solution does not depend on the g-value. 

As before, on-site surface solutions increase their width 
while the g- values increase [see Fig flTf inset]. For g < 0, 
the opposite is true. On the other hand, inter-site mode 
profiles are not really affected by the long-range term. In 
general, our analytical estimates for the large frequency 
limit of bulk solutions are the same than the ones for 
the surface solutions. In this limit, on-site solutions are 
composed essentially by one peak (independent of the 
particular excited site), while inter-site modes are com- 
posed by just two peaks. The SLM approximation for the 
on-site solution located at the right surface reduces to the 
dimer model with the same type of solutions previously 
described in section IV Therefore, a similar transition 
is expected for surface solutions; i.e., the on-site unstag- 
gered surface solution disappears for g > q while its fre- 
quency threshold also increases as g — » g, being infinite 
for q = g. For g < 0, the frequency /norm threshold for 
the on-site mode decreases as g becomes more negative. 



Fig|l3] shows a diagram with the numerically computed 
minimal norm (Nth) to excite a localized on-site solution 
at the surface for different values of g. We clearly see the 
effect of the long-range dipolar interaction: for g > the 
norm threshold increases, diverging for g — »• 1. This is 
because the first nearest-neighbor (U2) becomes more and 
more important and a localized solution requires now a 
larger norm (nonlinearity) to be trapped at the surface of 
the optical lattice. However, for g < the dipolar inter- 
action reduces the effective nonlinearity and, therefore, 
the frequency of the solution decreases, making possible 
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FIG. 13: Norm threshold versus dipolar interaction strength 
for q = 1. 

to excite a solution with smaller norm. These numerical 
results coincide perfectly with our estimates stemming 
from the dimer phenomenology. From the application 
point of view, a negative dipolar interaction would make 
possible the excitation and observation of surface states 
in simpler experimental conditions. Phenomenologically 
speaking, the effective energy potentials in the presence 
of a boundary look as the ones shown in Ref. [20] . Above 
norm thresholds, stationary solutions (extrema) can be 
excited in all sites of the lattice. However, for a smaller 
norm, the repulsive effect of the surface is non-negligible, 
and the nonlinearity is not able to create a localized state 
at the first or next sites of the lattice. 



VIII. CONCLUSIONS 

In this work we studied the interplay of scattering 
length and dipolar interactions on nonlinear localized 
modes (bulk and surface) of Bose-Einstein condensates 



in deep optical lattices. We started analyzing the pro- 
cess of modulational instability of nonlinear plane wave 
in a dipolar nonlinear lattice and established the regions 
of instability. In particular, for vanishing local atomic 
scattering length we showed that, for attractive dipolar 
interactions the nonlinear plane wave is linearly stable 
below a threshold amplitude, otherwise they are always 
unstable. That allowed us to find the favorable condi- 
tions for the existence of discrete solitons. Then, us- 
ing this information, the existence and stability of bulk 
discrete solitons was investigated analytically and con- 
firmed by numerical simulations. To study the existence 
and properties of surface discrete solitons we started with 
an analysis of the dimer configuration, where the prop- 
erties of symmetric and antisymmetric modes including 
the stability diagrams and bifurcations was investigated 
in closed form. For the case of a bulk medium, we an- 
alyzed properties of fundamental on-site and inter-site 
localized modes. In particular we showed that for attrac- 
tive local and nonlocal nonlinear it ies, there is a whole re- 
gion region in parameter space where both fundamental 
modes are simultaneously unstable. We also predicted 
a possible regime where the mode changes periodically 
from an inter-site mode (zero velocity), to an intermedi- 
ate state (maximum velocity) and then back to an on-site 
mode (zero velocity). This system could be used as a pre- 
cise 'atomic clock'. Results for expected shape of profiles 
agree well with the computed numerical profiles. We also 
investigated on-site and inter-site modes localized at the 
surface of ID lattice, i.e. surface solitons. We found 
and analyzed the form and properties of surface local- 
ized solutions finding a strong dependence on the effect 
of nonlocal nonlinear it ies. 
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